% This function computes the standard errors for the GMM estimates in
% the SEVAR(1) model when accounting for measurement errors in the factors
% and macros entering in the factor dynamics

function out = getAVarTheta2_threshold_macro(theta2,VarFactors,Cov10Factors,Cov01Factors,...
    x2,macros,bandwidthT,econAct,thresholdLevel,h0RegimeOn,hxRegimeOn)

% Computing the Jacobian matrix
epsValue    = 1D-7;
factors     = [macros;x2];
nx          = size(factors,1);
namesTheta2Full = fieldnames(theta2);
if h0RegimeOn == 1 && hxRegimeOn == 1
    % the full model
elseif h0RegimeOn == 0 && hxRegimeOn == 1
    % Remove h0_2
    for i=1:nx
        theta2 = rmfield(theta2,['h0',num2str(i),'_2']);
    end
elseif h0RegimeOn == 1 && hxRegimeOn == 0
    % Remove hx_2
    for i=1:nx
        for j=1:nx
            theta2 = rmfield(theta2,['hx',num2str(i),num2str(j),'_2']);
        end
    end
elseif h0RegimeOn == 0 && hxRegimeOn == 0
    % Remove h0_2
    for i=1:nx
        theta2 = rmfield(theta2,['h0',num2str(i),'_2']);
    end
    % Remove hx_2
    for i=1:nx
        for j=1:nx
            theta2 = rmfield(theta2,['hx',num2str(i),num2str(j),'_2']);
        end
    end
end
namesTheta2 = fieldnames(theta2);
numTheta2   = size(namesTheta2,1);
moments     = Var1Threshold_GMM_moments_macro(theta2,VarFactors,Cov10Factors,Cov01Factors,...
                x2,macros,econAct,thresholdLevel,h0RegimeOn,hxRegimeOn);
T = size(moments,2);
%Q = mean(moments,2)'*mean(moments,2)
R = zeros(numTheta2,numTheta2);
for i=1:numTheta2
    theta2Grad = theta2;
    theta2Grad.(namesTheta2{i}) = theta2Grad.(namesTheta2{i}) + epsValue;
    moments_p_eps = Var1Threshold_GMM_moments_macro(theta2Grad,VarFactors,Cov10Factors,Cov01Factors,...
        x2,macros,econAct,thresholdLevel,h0RegimeOn,hxRegimeOn);

    theta2Grad = theta2;
    theta2Grad.(namesTheta2{i}) = theta2Grad.(namesTheta2{i}) - epsValue;
    moments_m_eps = Var1Threshold_GMM_moments_macro(theta2Grad,VarFactors,Cov10Factors,Cov01Factors,...
        x2,macros,econAct,thresholdLevel,h0RegimeOn,hxRegimeOn);
    
    R(i,:) = (mean(moments_p_eps(:,1:T-1),2)-mean(moments_m_eps(:,1:T-1),2))'/(2*epsValue);
end

% The Newey-West estimate of S with lags controlled by bandwidthT (Robust to auto-correlation and
% heteroskedasticity)
GAMMA0 = zeros(numTheta2,numTheta2);
for t=1:T-1
    GAMMA0 = GAMMA0 + moments(:,t)*moments(:,t)'/(T-1);
end
S      = GAMMA0;
GAMMAk = zeros(numTheta2,numTheta2,bandwidthT);
for k=1:bandwidthT
    for t=1+k:T-1
        GAMMAk(:,:,k) = GAMMAk(:,:,k) + moments(:,t)*moments(:,t-k)';
    end    
    GAMMAk(:,:,k) = GAMMAk(:,:,k)/(T-1-k);
    S = S + (1-k/(1+bandwidthT))*(GAMMAk(:,:,k) + GAMMAk(:,:,k)');
end

% The asympotic covariance matrix 
invS             = S\eye(size(S,1)); %inv(S); 
tmp              = R*invS*R';
invTmp           = tmp\eye(size(tmp,1));
%out.VarRobust   = 1/(T-1)*inv(R*invS*R'); 
out.VarRobust    = 1/(T-1)*invTmp; 
out.namesFull    = namesTheta2Full;
out.names        = namesTheta2;
